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ABSTRACT 

In order to understand the influence of magnetic fields on the propagation 
properties of waves, as derived from different local helioseismology techniques, 
forward modeling of waves is required. Such calculations need a model in magne- 
tohydrostatic equilibrium as initial atmosphere to propagate oscillations through 
it. We provide a method to construct such a model in equilibrium for a wide range 
of parameters to be used for simulations of artificial helioseismologic data. The 
method combine the advantages of self-similar solutions and current-distributed 
models. A set of models is developed by numerical integration of magnetohydro- 
static equations from the sub-photospheric to chromospheric layers. 

Subject headings: MHD; Sun: magnetic fields; Sun: sunspots 



1. Introduction 



In the recent years, local helioseismology has provi ded new insights into the sub - photo 



spheric structure of quiet and active regions of the Sun (Duvall et al. 



1993 



Kosovichev 1999 



20021 : iKosovichev et all I2OO0I : Izhao fc Kosovichev! I2OO3I : braun fc Lindsevi boooh . However 



the influence of magnetic field s on data interpr e tation has not been fully e xplored. Theoreti- 



cal eff o rts have been made by Crouch fc Callv (hoOof) , 



f l2006f ). ICallv fc Goossensl (120071 ) 



Schunker et al 



Callvi fl2005l . l2006h . lSchunker & Callv 



J2008h m order to include mode conversion 



and to model ray path of waves in magnetized structures by means of analytical theory. These 
studies confirm the potential importance on helioseismic measurements of the so called sur- 
face effects caused by the presence of a magnetic field. A more complete understanding of the 
problem will, probably, be better reached via direct forward modeling of helioseismological 
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, e.g., 


Gizon et al. 


f 


2006); 


Khomenko & Collados 


f 


2006 


); 


Shelvag et al. 


f 


200' 




Hanasoge 


f 


2008 


); 


Cameron et al. 


i (except 


Shelvae et al. 


2007 


) , the authors apply a similar 



strategy. In particular, they assume the existence of an equilibrium atmosphere containing 
a magneto-static structure whose properties may resemble to a larger or lesser extent those 
of a sunspot or a magnetic flux tube. A small-amplitude perturbation is then applied to 
the system in order to study wave propagation, wave mode transformations, amplitude and 
phase behaviour of waves in a complex magnetic field topology, etc. 

The behaviour of waves observed in active regions is very sensitive to their magnetic 
field configuration. Photosphere and low chromosphere are regions where a small change of 
parameters such as the size of the magnetic structure, or its temperature, density or magnetic 
field strength and inclination, may produce significant changes in the resulted wave field. 
Simulations can be of invaluable help to explore, within a full parameter space, different 
magneto-static structures in order to understand the effects produced by the magnetic field 
on the measurable variables used in local helioseismology. 

The latter task requires a robust procedure to construct magneto-static structures of 
desired properties. In this paper, we propose a strategy with that aim and apply it to obtain 
thick structures, as prototypes for solar spots and pores. As a minimum requirement, the 
model should fulfill the following properties: (i) in the photosphere the model should, on 
average, reproduce the properties of a typical sunspot; (ii) at the border, the model should 
smoothly merge into a quiet-Sun non-magnetic model atmosphere; (iii) there should be a 
possibility to choose the profile of thermodynamic parameters at the sunspot axis; Wilson 
depression should be taken into account; (iv) magnetic field strength, inclination and the 
radius of the structure should be adjustable; (v) the model should be easily extensible into 
an arbitrary depth below the photosphere. 

There is a vast amount of works on magneto-static models reported in the literature. 
Leaving apart small-scale flux tube m odels, those for thick structures can be divided into 
those possessing a current sheet (e.g. iPizzol Il990f). with a sharp magnetic-non magnetic 



interface, and those with distributed currents (e. q. lPizzolll986l ). showing a smooth transition. 
Without discussing advantages and disadvantages of the both, we will proceed here with 
current-distributed models. 



Present current-distributed models apply two different philosophies. In the first cate- 
gory of models, the magnetic structure is prescribed and the distribution of thermodynamic 
variables is looked fo r to be in agreement with this structure. This is the class of self- s imilar 



models, proposed by ISchliiter fc Temesvaryj (119581 ) and then extended by, e.g., ILowi (11975 
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19801 ). In the second category, the pressure distribution is prescribed as the boundary con- 
dition at the axis of the magnetic structure and in the far-away non-magnetic atmosphere. 
Both, pressure and magn etic field, ar e iteratively changed in the remaining points to reach 
an equilibrium situation (|Pizzolll986l ). 



From the point of view of the requirem ents se t abov e, both classes of models have advan- 
tages and disadvantages. The approach by iPizzol (119861 ) is more fruitful in the photosphere, 
since the pressure distributions of the field-free and magnetized atmospheres can be taken 
from observations and are relatively well known. However, for deep sub-photospheric lay- 
ers, the models that can be taken as boundary conditions are scarce. More precisely, the 
quiet-Sun non-magnetic pressure stratification can be taken from helioseismologi cal data, 
for example, from the standard solar model of IChristensen-Dalsgaard et al . (1996) . As for 



the sunspot axis, no precis e data are availa ble (see, how ever I Zhao et al 



2002: Couvidat et al. 



2001 



Kosovichev 



20061 ). The model of IPizzol (119861 ) turns out to be very sensitive to 



the pressure deficit inside sunspots and the method, in general, has poor convergence if the 
simulation box is too deep and is very sensitive to the guess of the pressure distribution at 
the sunspot axis. This makes the method unsuitable for the purpose of our work. 



On the other side, the procedure proposed by iLowl (Il980l ) works better in deep layers, 
where the gas pressure dominates over the magnetic pressure. In the photosphere, where the 
plasma becom es magnetic ally-dominated, negative pressures are frequently obtained from 
the method of lLowl (Il980l ). It is complicated to guess the parameters of the magnetic field 
configuration in order to avoid this problem. At the same time, if one wishes to extend the 
models into the photosphere and higher layers, the magnetic field strength is l imited to rather 



low fl ux tube-like values, not appropriate for sunspots (see iHanasogd 120081 ; ICameron et al. 
2008h . 



In this paper, we take advantage of both Pizzo-like and Low-like approaches and propose 
a method to calculate the magneto-static equilibrium of a thick sunspot-like structure with 
the properties defined above. Below we describe the equations that allow to successfully 
merge results from both methods and show examples of MHS solutions for a wide range of 
parameters. The conclusions are given in the last section. 



2. Method 

We solve the equilibrium force balance equation together with divergence-free condition 
for the magnetic field: 

- VP +« + i(Vx5)xB = 0, (1) 
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VB =0. 



Following iPizzd (119861 ). the equations are solved in cylindrical coordinates (r, 0, z) and axial 
symmetry is assumed (i.e., all variables are independent of <fi). Under these conditions the 
magnetic field vector can be conveniently written in terms of the field line constant u: 

1 du G(u) 1 du 
r q z ' r ' r Q r 



B 



(2) 



where G(u) is a func tion related to t he twist component of the field. The va riable u is used 
both in IPizzd (119861 ) and iLowl (ll980l ). The di f ferenc e is that, in iLowl (ll980l ). the analytical 
expression for u is postulated, while, in IPizzd (119861 ). the shape of the field lines u is looked 
for by solving iteratively the equation of the force balance for given boundary conditions in 
agreem ent with som e pressure distribution. Except for constants, the functional form of u 
used by IPizzd ( 119861 ) as initial condition at th e low e r bou ndary of the computational domain 
is exactly the same as the one postulated by ILowl (jl980l ). Thus, both models can be joined 
in a natural way, assuming that the deep layers of the model sunspot can be approximated 
by the self-similar solution and the upper layers by the solution of Pizzo. 



In ILowl (ll980l ). following the spirit of self-similar solutions, the field line constant u is 
expressed SIS db function of one variable, <p: 



u(r, z) = u((p); (f = r-F(z); F(z) = (z 2 + a 2 ) 



2\-l 



(3) 



where a is a constant parameter. The field is untwisted and the azimuthal component, B^, 
is zero. This is equivalent to setting G{u) = in Eq. [21 Using the above expression, Eq. [2] 
can be rewritten as a function of p: 



B 



dz dcp' ' dip 



(4) 



Following Low, the function du/dip has to satisfy certain normalizations in order to fulfill 
the force balance equation. This leads to the following expression: 



= B^h 2 ■ exp(-rip) . 



(5) 



Here, Bq is a parameter that controls the magnetic field strength and h is a suitable length 
scale (note that, in the original paper of Low, dimensionless variables are used, while here 
we choose to use physical dimensions for all the variables). 

Introducing Eq. [5] into Eq. HI the horizontal and vertical components of the magnetic 
field vector in the Low's model are written as: 



B r (r, z) 



2#o L 



[z - z d )rh 2 



{{z - z A ) 2 + a 2 ) 2 



cxp 



(z - z d ) 2 + a 2 



(6) 
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B z (r, z) = 25 L — — exp ( , ^ ) (7) 



[z - z A ) 2 + a 2 \(z - z d ) 2 + a 2 / 
The parameter is a reference height, where the magnetic field is purely vertical. The latter 
expression is directly comparable with the one used in Pizzo as a boundary condition at the 
bottom boundary of the domain: 

,2- 



r 

e 



B z (r,z ) = B^ex P [- r -\ (8) 



Comparing these two expressions, we see that, if both models are to be joined at some 
arbitrary height, z = zq, the parameters of the models should be related as: 

p l 2/r 2 

Oo - z A ) 2 + a 2 

r 2 e = ((zo - z A f + a 2 )h (10) 
Keeping this in mind, the model can be constructed following the steps described below. 



2.1. Step 1: Generation of a self-similar solution in deep layers 



In deep sub-photospheric layers, we calculate a self-similar solution for B after Eqs. Eland 
[7J The pressure and density distributi ons with h eight and radius are found from analytical 
expressions (equations 50 and 51 in iLowl 119801 ) . As a boundary condition at the right 
boundary (field-free atmosp h ere) w e take the pressure and density from the model S of 



Christensen-Dalsgaard et al.l (119961 ). The choice of the field-free pressure and density is 



rather arbitrary, however. Because of azimuthal symmetry, the left boundary of the domain 
corresponds to r = 0, i.e., the axis of the magnetic structure. The lower boundary is taken 
exactly at height z A , where B r (r, z A ) = at all distances r. In all the cases presented here, 
Zd = — 10 Mm and the origin for the z-axis is taken at the base of the photosphere. Given 
pressure and density, the temperature distribution in the model sunspot can be calculated 
using the equation of state either in tabular form or the one for an ideal gas. The parameters 
7], a and Bq can be chosen freely. An additional free parameter is the height, zq, that 
limits the upper boundary of the self-similar model. Depending on the height of the upper 
boundary, Bq can be larger or smaller in order to prevent from getting negative gas pressures. 

The basic topology of the solution is given in Fig. [TJ for the parameters indicated in the 
figure caption. These parameters are chosen on purpose in order to demonstrate that the 
method is able to deal with lar ge field str engths. The sunspot radius at the bottom boundary 
is roughly defined by ari~ l l 2 ( Lowlll980l ). The inclination of the field at the top boundary 
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changes with the distance from the axis, from to about 70 degrees at the right-most point 
of the domain. The magnetic field is concentrated inside the first 15 Mm from the axis, 
being weak in the rest of the domain. The field strength drops at the axis from 12 kG at 
z = —10 Mm to 4 kG at z — — 1 Mm depth. The gas pressure is always above the magnetic 
pressure. The Wilson depression is rather weak and the model sunspot is almost thermally 
plane-parallel in deep layers, as follows from the distribution of the acoustic speed, cs- We 
do not have criteria to judge how realistic is this description. Data of the sub-photospheric 
distribution of the sound s peed in sunspots are scarce and uncertain (see, however , time- 
distance analysis results by lZhao et al.ll200ll ; iKosovichevi 120021 ; ICouvidat et al.ll2006l ). 



2.2. Step 2: Generation of potential solution in the overlaying atmosphere 

Given the values of r], a, Bq, z&, and zq we calculate the initial parameters of the Pizzo 
model from Eqs. M and [TUJ This giv es us Bp = 4 kG and r e = 9.4 Mm. We follow the 



same steps as in the original paper of iPizzol (119861 ) and start from computing the potential 
solution: 

d 2 u 1 du d 2 u 



dr 2 r dr ^ dz 2 

The bottom boundary of the domain coincides with the top boundary from the previous 
step and is located at z = — 1 Mm, below the photosphere. The field line constant u at the 
bottom boundary is approximated by: 

M = r e 2 <(i-exp(-^)/2 (12) 

At the left (sunspot axis) and top boundaries u = (vertical field) and u approaches a 
constant value at the right boundary (horizontal field). With this set of boundary conditions, 
the boundary value problem posed by Eq. [IT] is solved by standard methods. 



2.3. Step 3: Generation of magneto-static solution in the overlaying 

atmosphere 

The potential solution obtained in Step 2 is used as initial guess in the integration of 
the complete force balance equation along the magnetic field lines (equation 4 in the paper 
by Pizzo): 

d 2 u ldu^d 2 u ^ 2 dP(u,z) 
dr 2 r dr dz 2 du 
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In order to start iterations we need an app roximation fo r th e distributio n of pressure 
along the magnetic field lines P(u, z). Following iPizzol (119861 ) and iLowl (119751 ) we take it of 
the form: 

dz 



P(u, z) = Pq(u) exp 



(14) 



h(u, z) 

where Pq{u) is the gas pressure along the bottom boundary. The function h(u,z) is a scale 
height. For a complete description of the problem, the representative pressure distributions 
along the axis and in the field-free quiet atmosphere need to be s pecifie d. As field-free at- 
mosphere, we use of the model S of IChristensen-Dalseaard et al.l (119961 ). smoothly joint to 



the VAL-C m odel of the solar chromosphere (jVernazza et al.lll98ll ). At th e axis, we use the 



Avrettl (Il98ll ) model in the upper layers. In deep layers, we take a model by lKosovichev et al. 



(120001 ). obtained from helioseismic inversions of the sound speed beneath sunspots. This 
model already takes into account the Wilson depression, which is about 450 km. The com- 
plete model on the axis can be shifted up or down, though, if smaller or larger values of the 
Wilson depression are required. The model has a cool region just below the surface and a 
hot one below, down to about -10 Mm. For our purposes we take this model starting from 
-1 Mm depth and, thus, the hot layer is not taken into account. Once these models are 
specified, we calculate a smooth transition be tween them for the gas pressure P(u,z) and 
scale height h(u,z) distributions, as given by IPizzol (119861 ) in his equations (13), (17) and 
(18). Then, the Eq. [13] is iterated until a convergence criterion is reached. 

The model sunspot for the parameters given above is shown in Fig. [2j The magnetic 
field on the axis drops from 4 kG at z — — 1 Mm to about 2 kG at = 2 Mm, which is a rather 
large value at this height. Due to such a large field strength, the Alfven speed exceeds 10 4 
km s -1 in the upper layers. The image of the sound speed shows the presence of the Wilson 
depression around z = 0, i.e. the temperature at the sunspot axis is smaller than in the 
outside atmosphere at a given height. Note that at higher layers, the effect is the opposite 
and the temperature inside the sunspot is larger. This effect is due to initial distribution 
in the model atmospheres taken as boundary conditions. The field lines are more inclined 
comparing to the Low solution in Fig. [TJ 



2.4. Step 4: Concatenating the solutions 

Both solutions obtained in Step 1 and Step 3 are in MHS equilibrium. In order to 
construct the complete model from deep to high layers one has to put one model on top of 
the other. However, despite B z at the bottom boundary of the Pizzo model is calculated in 
agreement with B z at the top boundary of the Low model, there is a discontinuity in the 
horizontal component of the magnetic field. This discontinuity can be appreciated in Fig. [31 
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where the field lines with a dashed style correspond to the two models concatenated as they 
are. The reason of this discontinuity is two-fold. On the one hand, the physics of the solution 
changes abruptly from one model to another, thus changing the gradients of the magnetic 
field, gas pressure, etc. On the other hand, the boundary condition for the field line constant 
u is not the same in the both models. In the case of Low model, there is no need to put a 
condition on u, neither there is a possibility. The inclination of the magnetic field lines at the 
right boundary is a consequence of the parameters of the model and should not be necessarily 
horizontal. Contrarily, in the case of the Pizzo model, we impose horizontal magnetic field at 
the right boundary. The dependence of B z on r given by Eq. [8] defines the vertical magnetic 
field strength but does not put constraints on the horizontal field component. 

Thus, in order to obtain a smooth solution everywhere in the domain, we repeat the 
Step 3 calculations for the complete model sunspot. We take the pressure distributions at 
the axis and in the field-free outside atmosphere from the joint model at all heights. The 
boundary conditions for u are the same as in the Pizzo model. The distribution of u at the 
bottom boundary is taken from the Low model. Then we repeat the solution of Eq. [13j The 
resulted topology of the magnetic field lines is plotted in Fig. [3] in solid line style. The field 
lines in the final solution are more horizontal in low layers, while they are more vertical in 
the upper layers. 

Fig. H] gives the distribution of some parameters in the complete model sunspot, at all 
layers. We can see that the last iteration has re-distributed all the parameters compared to 
the individual Low and Pizzo parts of the solution. In particular, the magnetic field gradient 
at the axis is now more steep and the field strength at high layers becomes lower. The field is 
in general more inclined, being horizontal at the right hand side domain boundary, consistent 
with our imposed boundary condition there. The gas pressure is modified accordingly to 
maintain the new force balance. 

Fig. [5] gives a more detailed view on the model spot solution. It shows the distribution 
with radius (left panels) and with depth (right panels) of some parameters of the sunspot 
atmosphere. The field strength decreases rapidly with height at an average rate of about 
1 G/km at the axis. The magnitude of the gradient decreases with height and with dis- 
tance to the axis. These g radients are in agreement with photospheric spectropolarimetric 



observations (jSolankill2003l ). The magnitude of the pressure deficit inside the model sunspot 
decreases with depth almost disappearing at about -2 Mm depth, in accordance with our 
assumption of self-similarity of the MHS solution at larger depths. 

As can be seen from the radial pressure distribution, there is a pressure excess observed 
at larger heights in the chromosphere at some distance from the axis. This pressure excess 
would produce a bright ring in the emergent intensity from the model sunspot and it is 



- 9- 



present in the original model by iPizzol (119861 ). As shown in the latter work, the bright ring 



can be removed by an improved initial guess of the P(u, z) distribution. It is unimportant 
for the purpose of the present work since we only need an approximate agreement between 
the average properties of the MHS solution and the observed properties of sunspots. The 
magnetic field lines of the model sunspot are inclined less than 30 degrees within the first 10 
Mm from the axis, which can be considered as the umbra. Due to the boundary condition, 
the inclination changes gradually becoming 90 degrees at the edge of the model, where the 
magnetic field is already very weak. The ratio between the sound speed and the Alfven 
speed squared (which gives the measure of the gas to magnetic pressure) changes orders of 
magnitude from 10 6 at z — —10 Mm to 10~ 6 at z =2 Mm. Note that despite this, there is 
no problem with the convergence of the solution. 

In the next section, we give more examples of MHS solutions comparing models obtained 
with various sets of parameters. In the examples below we discuss the models calculated in 
a complete domain from z = — 10 to z = 2 Mm. 



3. Examples 

Dependence on magnetic field strength. Fig. [H] shows the magnetic field topology 
of models with different values of the magnetic field strength (parameter Bfr), all the other 
parameters being exactly the same as previously. As can be seen by comparing the different 
curves on the figures, the resulting gas pressure stratifications only differ in the highest layers 
by the amount of the pressure deficit. The magnetic field topology is indistinguishable in 
all the cases. This property originates from two effects. On the one side, the self-similar 
solution in the bottom part of the domain scales with magnetic field strength, i.e., the field 
line topology does not depend on Bq . On the other side, the Pizzo solution in the upper part 
is close to potential, imposed by the solution at the bottom part. The potential solution also 
scales with the magnetic field strength and is independent of the thermodynamic properties. 
Both effects lead to the independence of Bq of the magnetic topology of the final solution 
in the complete domain. This is a useful property from the point of view of helioseismology 
simulations. Using a set of models with different magnetic field strength, but otherwise the 
same, the effects of the magnetic field strength on waves can be checked independently of 
the effects of the magnetic field inclination. 

It should be noted that the above property originates only from the particular choice of 
the parameters a and rj. This choice produces r e large enough, so that the final solution in 
the upper part of the domain ap proaches to potential and becomes almost independent on 



the pressure distribution P(u, z) (IPizzol 1 1986h 



-10- 



The models presented in this Section are available electronically in FITS format from 
the Astrophysical Journal web site. 

Dependence on a and 77. Fig. [7J shows the magnetic field topology of the models 
calculated with different values of the parameter 77 (Eqs. [6] and [7]), the rest of the parameters 
being the same. Note that, according to the above equations, the inclination of the magnetic 
field is independent of rj in the Low model at the bottom part of the domain. However, the 
initial radius of the structure in the Pizzo part of the solution r e (Eq. [TUT) depends on 77, 
thus changing the inclination of the magnetic field lines in the upper part of the atmosphere. 
The final iteration performed in Step 4 takes that into account, making the solution in the 
complete domain dependent on 77. 

The change of 77 produces two effects. By increasing 77, we decrease the magnetic field 
strength by a smaller amount than by varying Bq, as in the previous example. At the same 
time increasing 77 produces an increase of the inclination of the magnetic field lines in the 
solution in the complete domain. The difference in the inclination is more pronounced in 
the deep layers of the model. The magnetic field topology of the solutions is different. The 
gradient of the magnetic field at the axis is slightly larger for larger values of 77 in the sub- 
photospheric part of the model. Again, this difference is produced after the final iteration 
in Step 4 since the B z (r = 0) given by Eq. [JJof the Low part of the solution is independent 
of 77. The pressure distribution at the axis and the amount of the pressure deficit are not 
very different between the given models. The models presented in this Section are available 
electronically in FITS format from the Astrophysical Journal web site. 

Varying the parameter a produces similar effects. The difference is that, by varying a, 
we change mostly the curvature of the magnetic field lines and the radius of the structure, 
not affecting much the magnetic field strength. 

Dependence on zq. Another parameter introduced in our modeling is the height 
where both solutions merge, zq. Fig. [8] shows the topology of the magnetic field lines and 
the pressure and density distributions with height at the axis, for three models with different 
values of z . In this example we take different values of a and 77 compared to above examples 
(see figure caption) in order to produce a structure with a smaller radius. This way, we 
show that the procedure is robust and can produce magnetic structures with very different 
properties. 

The magnetic field strength at the axis is almost independent of the choice of Zq. The 
amount of the pressure deficit at the near-surface layers increases with decreasing z from 
— 1 to —3 Mm, extending to larger depths. Note, however, that we can not shift the level 
of zq much deeper than —3 Mm due to a poor convergence of the solution. Despite the 
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magnetic field strength is nearly the same, the position of the (3 — 1 level is different in all 
the solutions due to the different amount of the pressure deficit. 

The inclination of the magnetic field lines is similar in the central part of the model 
sunspots in the three solutions. At the periphery, especially at larger depths, the field lines 
get more inclined with decreasing zq. Thus, the field is more concentrated to the central part 
and the effective radius of the structure is smaller. The models presented in this Section are 
available electronically in FITS format from the Astrophysical Journal web site. 

4. Conclusions 

In this study, we propose a method to construct a magnetostatic structure with prop- 
erties and size of a typical sunspot, from the deep interior to the solar surface. Previously 
published methods to construct such a model fail due to a poor knowledge of the thermo- 
dynamic and magnetic parameters of sunspots in sub-photospheric layers. We make use 
of self-similar models in deep layers and show that such models can naturally merge with 
models where the pressure distribution is prescribed on the well as the field-free 

atmosphere, allowing for a more realistic description of the atmospheric layers of sunspots. 
The procedure shows a rather good convergence and stability. By changing the parame- 
ters of the solution, a set of models can be produced with desired properties. We suggest 
that these models may be used, among others, in artificial helioseismology data simulations. 
Given a set of models, a parametric study can be done, investigating the influence of the 
topology and strength of the magnetic field of sunspots on the parameters inferred by the 
local helioseismology measurements in solar active regions. 

This research has been funded by the Spanish Ministerio de Education y Ciencia through 
projects AYA2007-63881 and AYA2007-66502. 
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Fig. 1.— Topology of the Low solution with a = 2h, h = 3 Mm, r? = 1.3, B = 25000 G 
and Zq = —1 Mm. Top: magnetic field strength; middle: acoustic speed; bottom: log of the 
Alfven speed. White lines are magnetic field lines. Black lines with labels are the contours 
of the ratio of the sound speed and the Alfven speed squared, c 2 s /v\. 
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Fig. 2. — Topology of the Pizzo solution with Bq = 4 kG and r e = 9.4 Mm. Top: magnetic 
field strength; middle: acoustic speed; bottom: log of the Alfven speed. White lines are 
magnetic field lines. White lines with labels are the contours of Cg/v\. Note that, for better 
visualization, the vertical axis has been expanded. 
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Fig. 3. — Topology of the magnetic field lines before iterations (dashed) and after a new 
equilibrium is reached (solid). 
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Fig. 4. — Topology of the complete solution with Bq = 25000, a = 2h, h = 3 Mm and 
r) = 1.3 (Bq = 4 kG and r e = 9.4 Mm). Top: magnetic field strength; middle: acoustic 
speed; bottom: log of the Alfven speed. White lines are magnetic field lines. White lines 
with labels are the contours of c$/v\. 
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Fig. 5. — Distribution with radial distance (left panels) and with depth (right panels) of 
the magnetic field strength, pressure, ratio c 2 s /v 2 A and the magnetic field inclination for the 
sunspot with Bq = 25000, a = 2h, h = 3 Mm and rj — 1.3 (corresponding to Bq = 4 kG 
and r e = 9.4 Mm at z — — 1 Mm). The radial pressure distributions are normalized to their 
values at the right boundary. 
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Fig. 6. — Top panels: height dependence of the magnetic field and pressure at the axis for 
the models with a = 2h, h = 3 Mm, r? = 1.3 and B^= 10000 (red line), 16000 ( green line), 
25000 (blue line) and 40000 (magenta line) G. Bottom panel: topology of the magnetic field 
lines for the same solutions (same color coding). Contours of the magnetic field strength of 
B = 1000 G are shown by dotted lines for each case. Horizontal solid lines mark the levels 
of C S = Va- 
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Fig. 7. — Top panels: height dependence of the magnetic field and pressure at the axis for 
the models with a = 2h, h = 3 Mm, Bq= 25000 and rj = 1.3 (red line), 2.5 (green line) and 
3.5 (blue line). Bottom panel: topology of the magnetic field lines for the same solutions 
(same color coding). Contours of the magnetic field strength of B = 1000 G are shown by 
dotted lines for each case. Horizontal solid lines mark the levels of c$ = va- 
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Fig. 8. — Top panels: height dependence of the magnetic field and pressure at the axis for 
the models with a = 1.5h, h = 3 Mm, Bq = 25000, rj = 3.5 and zq=-1 Mm (red line), -2 Mm 
(green line) and -3 Mm (blue line). Bottom panel: topology of the magnetic field lines for 
the same models (same color coding). Contours of the magnetic field strength of B = 1000 
G are shown by dotted lines for each case. Horizontal solid lines mark the levels of cs = va- 



